//
// labFisStefBoltElab.cc
//
// codice per elaborazione dati esperimento Stefan Boltzmann
//
// Istruzioni
//
// il programma lavora su un file di dati "labFisStefBoltSave.root"
// per crearlo (se non c'e gia') occorre compilare ed eseguire
// "labFisStefBolt", poi cambiare il nome del file prodotto:
// # mv labFisStefBolt.root labFisStefBoltSave.root
//
// il programma produce un fine con i dati elaborati di nome:
// "labFisStefBoltSaveR.root".
//
// - per eseguirlo occorre far partire root:
// # root <return>
//
// - eseguire il programma da dentro Root
// root [0] .x labFisStefBoltElab.cc
//

#include "TPad.h"
#include "TFile.h"
#include "TTree.h"
#include "TProfile.h"

void labFisStefBoltElab()
{
  TFile* f = new TFile("labFisStefBoltSave.root");
  TTree* nt = (TTree *)f->Get("nt");
  Long64_t nentries = nt->GetEntries();

  //
  // i dati corrispondono a 1000 campionamenti per 9 voltaggi (0V..8V)
  // creiamo istogrammi dove ogni "bin" corrisponde a 1000 campioni
  // cosi' avremo in ogni "bin" la media delle quantita' istogrammate
  // (voltaggio su bulbo e su resistenza)
  //

  Int_t channels = static_cast<Int_t>(nentries / 1000);

  TProfile* hvb  = new TProfile("hvb", "tensione sulla lampadina", channels, 0, nentries);
  TProfile* hvs  = new TProfile("hvs", "corrente al variare del voltaggio", channels, 0, nentries);

  //--- fill histogram
  nt->Project("hvb", "vb:Entry$");
  nt->Project("hvs", "vs:Entry$");

  // hvb->DrawCopy();
  hvs->DrawCopy();

  gPad->Update();

  //--- get histogram parameters
  UInt_t bins = hvb->GetNbinsX();
  // Double_t xlow = hvb->GetXaxis()->GetXmin();
  // Double_t xupp = hvb->GetXaxis()->GetXmax();

  TFile* hfile = new TFile("labFisStefBoltSaveR.root",
			   "RECREATE", "ntupla mediata");
  
  //--- allocate Root Tree
  TTree* tree = new TTree("ntr", "Legge di Stefan-Boltzmann");
  
  Float_t ntVbSet(0);
  Float_t ntVb(0);
  Float_t ntVs(0);
  Float_t ntVbErr(0);
  Float_t ntVsErr(0);
  Float_t ntT1(0);
  Float_t ntT2(0);
  Float_t temp(0); //temperatura calcolata col polinomio
  // Float_t tempErr(0);
  Float_t ress(0); //resistenza calcolata col polinomio
  // Float_t ressErr(0);
 
  //--- add branches, each with a single leaf
  tree->Branch("vbSet", &ntVbSet, "vbSet/F");
  tree->Branch("vb", &ntVb, "vb/F");
  tree->Branch("vs", &ntVs, "vs/F");
  tree->Branch("vbErr", &ntVbErr, "vbErr/F");
  tree->Branch("vsErr", &ntVsErr, "vsErr/F");
  tree->Branch("t1", &ntT1, "t1/F");
  tree->Branch("t2", &ntT2, "t2/F");
  tree->Branch("temp", &temp, "temp/F");
  tree->Branch("ress", &ress, "ress/F");

  for(UInt_t i=1; i< bins+1; i++) {
    ntVb = hvb->GetBinContent(i);
    ntVbErr = hvb->GetBinError(i);
    ntVs = hvs->GetBinContent(i);
    ntVsErr = hvs->GetBinError(i);
    ress = ntVb/ntVs;
    temp=ress/1.5315;
    temp=103.898+214.93*temp-2.9944*pow(temp,2)+0.04328*pow(temp,3);

    // cout << ntVb << endl;
    tree->Fill();
  }

  tree->Write();
  hfile->Close();
  f->Close();
}
